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The quantum statistics (QS) correlations of identical bosons are well known to be sensitive to the 
space-time extent and dynamics of the particle emitting source in high-energy collisions. While two- 
pion correlations are most often experimentally measured, the QS correlations of three-pions and 
higher are rarely explored. A set of techniques to isolate and analyze 3- and 4-pion QS correlations 
is presented. In particular, the technique of built correlation functions allows one to more easily 
study the effects of quantum coherence at finite relative momenta instead of at the unmeasured 
intercept of correlation functions. 


I. INTRODUCTION 

In high-energy collisions, the quantum statistics (QS) 
correlations of identical particles provide the most sen¬ 
sitive probe of the ^ace-time structure of the particle- 
emitting source [ll-Q- Most often, 2-pion Bose-Einstein 
correlations are experimentally measured for which the 
techniques are well established. However, the statistical 
precision of modern experiments has opened the door to 
measurements of higher order QS correlations. 

Higher order QS correlations contain additional infor¬ 
mation about the source which cannot be learned from 
2-particle correlations alone. In particular, the suppres¬ 
sion expected from quantum coherence of identical pions 
in the final state increases substantially for higher orders. 
The possibility of quantum coherence in high en ergy col¬ 
lisions has been considered several times before (^-l9| . 

One such mechanism is the disoriented chiral conden¬ 
sate (DCC) The chiral condensate can be char¬ 

acterized with four scalar fields (cr, 7 ?) and is non-zero 
at low temperatures where chiral symmetry is sponta¬ 
neously broken. In the core of the collision, the chiral 
condensate is expected to vanish and the symmetry is 
restored. A hot expanding shell of particles essentially 
shields the core from the true vacuum outside the shell. 
As the core energy density drops, the chiral condensate 
forms and may be temporarily disoriented in one of the 
7? directions instead of the true vacuum cr direction. The 
disoriented vacuum eventually reorients in the cr direc¬ 
tion and is accompanied by coherent pion radiation. The 
DCC is very different from a Bose-Einstein condensate 
for which the latter develops from critically large pion 
densities. 

A quantum coherent or pure state of matter can be 
a rather delicate state of matter. Before the collision of 
two heavy-ions, each ion is in a separate pure quantum 
state. After the collision, the detected particles may be 
in a mixed state or even remain in a pure state depending 
on the number of particles detected (quantum entangle¬ 
ment) [TTj |. During intermediate times when the particle 
density is high, scatterings may also play a role of quan¬ 
tum decoherence. The hot and dense medium created is 
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often modeled with hydrodynamics which assumes local 
thermal equilibrium and quantum decoherence above the 
thermal length scale. Civen the success of hydrodynamic 
modeling of high-energy collisions, it is conceptually dif¬ 
ficult to imagine the survival of coherence in the final 
state. 

The effect of coherence manifests itself most direc^ 
in the suppression of Bose-Einstein correlations Mi- 
However, even for a coherent fraction as large as 50%, the 
intercept (vanishing relative momentum) of 2-pion corre¬ 
lations is decreased only by 25%. Furthermore, Coulomb 
repulsion and the dilution from long-lived resonance de¬ 
cays also suppress Bose-Einstein correlations. Civen the 
various sources of suppression and the unknown func¬ 
tional form of Bose-Einstein correlations, it is practically 
impossible to determine the coherent fraction from 2-pion 
correlations alone. 

In this article, a set of techniques for the measurement 
of 3- and 4-pion QS correlations is presented. In particu¬ 
lar, it is shown how the comparison of 2-, 3-, and 4-pion 
QS correlations can allow for a less ambiguous determi¬ 
nation of the coherent fraction. To that purpose, the 
method of constructed or built correlation functions is 
introduced. The techniques are presented for charged pi¬ 
ons up to 4*^ order but can be easily generalized beyond 
this. 

This article is organized into 10 sections. In Sec. H, 
the correlation functions and projection variables are de¬ 
fined. In Sec. HI, the concept of symmetrization and pair 
exchange amplitudes is introduced. The main technique 
of built correlation functions is introduced in Sec. IV. A 
short discussion of radii measurements using multi-pion 
correlations is given in Sec. V. Final-state-interactions 
are discussed in Sec. VI. Dilution effects from long-lived 
emitters are discussed in Sec. VH. In Sec. VHI, the cu- 
mulant and partial cumulant correlation functions are 
defined. A variety of model calculations are shown in 
Sec. IX. In Sec. X, the effect of multi-boson distortions on 
the comparison of built and measured correlation func¬ 
tions is calculated. Finally, the main findings of this anal¬ 
ysis are presented in Sec. XL 
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II. DISTRIBUTIONS AND CORRELATION 
FUNCTIONS 


The n-particle inclusive momentum spectrum is given 
by Nn{pi,p 2 , ■■■,Pn) where pi denotes the momentum of 
particle i. The product of n single-particle inclusive spec¬ 
tra, Ni{pi)Ni{p 2 )...Ni{pn), forms a combinatorial refer¬ 
ence sample to the full n-particle spectrum. Correlation 
functions of order n can then be constructed in the usual 
manner as the ratio of the two distributions, 


Cn{Pl,P2,---,Pn) 


Nn{pi,P2,---,Pn) 

Nl{pi)Ni{p2)...Ni{pn)' 


( 1 ) 


Experimentally, the two distributions can be measured 
by taking all n particles from the same collision event 
or by taking all n from different but similar events. In 
quantum optics, the maximum of Cn for identical pions 
is n! in the dilute gas limit [l^ - [l^ and neglecting isospin 

i- 


A. Projection variables 


Two-particle correlations are projected onto relative 
momentum in the Pair-Rest-Frame (PRF) or in the 
Longitudinally Co-Moving System (LCMS) where the 
z-component of the pair momentum vanishes. In the 
PRF, the ID Lorentz invariant relative momentum, q = 
\/—{Pi — P 2 )^{pi — P 2 )fi, is used. In the LCMS frame, 
the relative momentum vector is projected onto 3 dimen¬ 
sions. The projection onto the pair momentum vector 
forms (/out- The projection onto the longitudinal direc¬ 
tion (beam-axis) forms (/long- The direction perpendicu¬ 
lar to both “out” and “long” forms the gside projection. 

For 3- and 4-particle correlations, the Lorentz invariant 
relative momentum is used and defined by the quadrature 
sum of pair invariant relative momenta 


Qs — \j^12 T 9i3 T 923! 

<34 = \/ qh + 9i3 + 9i4 + Q23 + 924 + ^34 


( 2 ) 

(3) 


for 3- and 4-particle correlations, respectively. The trans¬ 
verse momentum dependence can be studied by also pro¬ 
jecting onto the average transverse momenta 


k'T = 


Kt,3 = 


At,4 = 


|pt,i + PT. 2 I 
2 ’ 

|Pt,1 -I- PT,2 + Pt.sI 
3 ’ 

|PT,1 +/?T.2 +PT.3 +PT. 4 I 


(4) 

(5) 

( 6 ) 


for 2-, 3-, and 4-particle correlations, respectively. 


III. SYMMETRIZATION 

For identical pions, the effect of quantum statistics is 
represented as a symmetrization of pion creation points. 


The symmetrization is only valid for independently pro¬ 
duced pions As pions in a coherent state are collec¬ 

tively produced, the symmetrization is absent for pairs 
of coherent pions. It is, however, present for mixed-pairs 
with one pion from the coherent pool and the other from 
chaotic pool of particles. 

In the Wigner function formalism and neglecting par¬ 
ticle interactions, the correlation functions are written in 
terms of the Fourier transformed single-particle emission 
function, S{x,p), which describes the phase-space distri¬ 
bution of particle production m. For the case of chaotic 
-|- coherent emission, the emission function is split into 
the sum of a chaotic and coherent parts BiQ 

S{x,p) = [I - G(p)]5'ch(a;,4') + G{p)Scoh{x,p). (7) 

The fraction of pions from the coherent pool is given by 
G{p) and is momentum dependent in general although it 
is treated as momentum independent for the rest of this 
article. The normalized pair exchange amplitudes, dy, 
are introduced in order to write the correlation functions 
compactly in terms of the emission functions 

f d^xSix,Kij)e'‘>ii^ 

[J <i*xS{x,pi) J d*yS{y,pj)]^/^ 

^ f d^x [(1 - G)Seh(x, Kjj) + GScoh(^, Kjj)] 

[f d'^xS{x,pi) f d‘iyS{y,pj)]^/^ 

= (1 - . ( 8 ) 

The magnitudes of the Fourier transform (FT), Ty and 
ty, are referred to as the pair exchange magnitudes for 
the chaotic and coherent parts of the source, respectively. 
The phases of the FT for the chaotic and coherent parts 
are given by <I>y and 4>ij, respectively. The average and 
relative pair momentum are given by Ay = {pi -\-pj)/2 
and Qij = Pi — Pj , respectively. The two-^ion correlation 
function can then be written as in Ref. 

C'2(pi,P2) = l (9) 

, \ f d^xS(x,Ki2)<^'n2^\'^-G^\ f d^xS^^i,(x,Ki2)e'n2^\^ 

” / d*xS(x,pi) f d*yS(y,P 2 ) ' ' 

= 1 -t (1 - -H 2G(1 - G)Ti 2 ti 2 cos(<I.i 2 - <> 12 ) (11) 

From Eq. [TI] one notices that while the FT phases 
disappear in the case of fully chaotic emission (G = 0), 
they are present in the case of partial coherence. Each 
term in Eq. [TT] can be represented in a symmetrization 
diagram similar to Ref. [I9l| and is shown in Fig.[TJ Pions 
from the chaotic pool are represented with a solid circle 
while those from the coherent pool are given by a hol¬ 
low circle. Chaotic pions yield a multiplicative factor of 
1 — G while coherent pions yield a factor of G. Solid and 
dotted lines yield a multiplicative factor from the pair ex¬ 
change amplitude for the chaotic (Ty-e***^ ) and coherent 
(tye*"^*^ ) part of the source, respectively. Lines running 
from pion j to pion i are the complex conjugate of those 
in the opposite direction. The sum of all graphs yield 
the full QS correlation function. Also, from Eq. |8] one 
has the equality, <I>y = and (j)ij = 

The symmetrization diagrams for 3- and 4-pion corre¬ 
lations are shown in Fig. [2] and [3l respectively. The 
number of permutations in the n-pion set is n!. From 
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• = 1-G (chaotic pool) 


o - G (coherent pool) 
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FIG. 1. 2-pion symmetrization diagrams. Derived from simi¬ 
lar figures in Ref. [l^ . 



C3 — 1 + 






+ 2 permutations 



+ 1 permutation 


FIG. 2. 3-pion symmetrization diagrams. Derived from simi¬ 
lar figures in Ref. [l^ . 


the diagrams in Fig. [JJ the S-pion correlation function 
can then be written as in Ref. Q 


In the absence of pion coherence (G = 0), the 3-pion 
FT phase factor in Eq. [TH (cos($i 2 + ‘^23 + ‘hai)) may 
be isolated by comparing 3-pion cumulant to 2-pion QS 
correlations. The 3-pion QS cumulant correlation has, 
by defini tion , all 2 -pion symmetrization terms removed 
from Eq. [T^] and is therefore given by 

03 ( 1 , 2, 3) = 1 -I- 2 T 12 T 23 T 31 cos('I>i 2 -I- ‘1>23 + 'I’ 3 l), (13) 

where the simplified notation (i) stands for pi. One may 
isolate the 3-pion phase factor by normalizing C 3 w ith 
lie app ropriate 2-pion correlation factors given in Eq. 1111 

illlai: 


r3(l,2,3) 


_C3(1,2,3) - 1_ 

V(C2(1,2) - 1)(C2(2,3) - 1)(C2(3,1) - 1) 

2cos(<l>i 2 -I->I‘23 +'I’ 3 l). (14) 


In Ref. [ 1 ^ it was shown that to leading order, the rel¬ 
ative momentum dependence of r^, is quartic. In the ab¬ 
sence of the 3-pion phase, is equal to 2.0 for all triplet 
relative momenta. The function was recently mea¬ 
sured at the LHC [2l| in Pb-Pb collisions. Although the 
systematic uncertainties were quite large for high Qa, no 
significant Q 3 dependence was observed for both low and 
high transverse momentum. 

Similarly, the 4-pion FT phase factor (cos(<f’i 2 + $23 + 
4*34 -I- (& 41 )) may be isolated by comparing 4-pion cumu¬ 
lant to 2 -pion QS correlations. The 4-pion QS cumulant 
correlation has all 2-pion, 2-pair, and 3-pion symmetriza¬ 
tion terms explicitly removed and is therefore given by 


<^3(pi, P2, P3) = 1 + (1 - G)2 (Tfa + TI 3 -h T|i) 

-I- 2G(1 — G)[ri2tl2 COs(<I>i2 — 012) 

+ 723423 COs($23 — 023) 

-t 73 i 431 COs(<I>3l — 03l)] 

-I- 2G(1 — G)^[Ti2T23431 COs('3?i2 -I- <I>23 + 03l) 

+ Tl2t23T3l COs(4>12 -I- 023 + 'I’3l) 

+ tl2T23T3l COs(0i2 + <I>23 + 'I’3l)] 

-I- 2(1 — G)®Ti2T23T3i COs('I>i2 + '3?23 + <i?3i).(12) 


04 ( 1 , 2, 3, 4) = 1 -I- 2 T 12 T 23 T 34 T 41 COs(<I>i 2 -I- '1>23 + 't ’34 + 'l’ 4 l) 

+ 2Ti2T24r43r3i COs(<i?i 2 -|- 'J>24 + ‘I ’43 + 'J>3l) 

+ 2Ti3r32r24r4i COs(<I>i 3 -I- 'J>32 + 'I>24 + 'J> 4 l) 

^1-1- 6 T 12 T 23 T 34 T 41 COs(<I>i 2 + <l>23 + ' 3?34 -I- <l> 4 i), (15) 

where the last equality is naturally obtained after aver¬ 
aging over many pion quadruplets. The 4-pion FT phase 
can then be isolated with r 4 defined as 

, , C4(1,2,3,4) - 1 

r4 1,2,3,4 = , y ^ ^ ^ / 

y(G2(l,2) - 1)(G2(2,3) - 1)(G2(3,4) - 1)(G2(4,1) - 1) 

= 6cos(4?i 2 -I-<l>23 +'1'34 +'I’4i). (16) 


C, = 1 + (^+(^)+0 

■66<’0^00^00<'0 

2 4 

:A:A"A"A 


+ 5 permutations 

+ 2 permutations 
+ 7 permutations 

+ 5 permutations 


FIG. 3. 4-pion symmetrization diagrams. Derived from simi¬ 
lar figures in Ref. [l^ . 


IV. BUILDING MULTI-PION QS 
CORRELATION FUNCTIONS 

In the absence of multi-pion FT phases and coherence, 
higher order correlation functions (n > 2 ) do not con¬ 
tain any additional information beyond that which is al¬ 
ready present in 2-pion correlation functions. The ab¬ 
sence/presence of both phenomena may then be tested 
by comparing the measured multi-pion QS correlations to 
the expectations from measured 2-pion correlations. For 
the rest of this section the FT phases (4’^ = 0, cj)ij = 0) 
are neglected as one expects a rather weak relative mo¬ 
mentum dependence |18j and since the recent ALICE 
results show no substantial Q 3 dependence [2l|. With 
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this simplification the building blocks of higher order 
correlation functions are the pair exchange magnitudes, 
and the coherent fraction, G. 

From Figs. [T][31 2-, i, and 4-pion QS correlations can 
be written as in Ref. [j: 

C2(l, 2) = 1 + (1 - G)‘^T^2 + 2G(1 - G)Ti2ti2 (17) 

C3(l, 2, 3) = 1 + (1 - Gf{Th + r|3 + r|i) 

+ 2G(1 — G)[Ti2tl2 + 723423 + 731 * 31 ] 

+ 2G(1 — G)^[Tl2T23t3l +Ti2*23731 + *12723T3i] 

+ 2(1 — G)®Ti2T23r3i (18) 

G4(1, 2, 3, 4) = 1 + (1 - Gf{Tl^ + r|3 + Til + 

+ 2G(1 — G)[Ti2*12 + 723*23 + 731*31 
+ ‘14’ + ‘24’ + ‘34’] 

+ (1 - Gf[Tl^Ti^ + T|iT|4 + T 24 TI 3 ] 

+ G(1 - G)3[Ti22T34*34 + TiiTl2tl2 
+ ‘13,24’ + ‘14,23’] 

+ 2(1 — G)^[ri2T23T3l + T 12 T 24714 
+ 731 T 34 T 41 + T 23 T 34 T 42 ] 

+ 2G(1 — G)^[Ti2T23*3l + Ti2*23 731 + *12 723731 
+ ‘124’ + ‘134’ + ‘234’] 

+ 2(1 - G)'‘[Ti2T23r34r4l 
+ 7’l2 724 743 r3i + Ti3T32T24r4i] 

+ 2G(1 — G)®[Ti2T23T34*i4 + ri2T23*34ri4 
+ 7 ’i2*237347’i 4 + *127237347l4 
+ ‘4321’ + ‘1243’ + ‘3421’ + ‘1324’ + ‘4231’], (19) 


where ‘ij\ Hj,kl\ ‘ijk\ and ‘ijkV represent permuta¬ 
tions within the same symmetrization sequence. In the 
case of fully chaotic emission, T^- can be extracted di¬ 
rectly from 2-pion correlations. As T^- is a 6D func¬ 
tion, experimental 2-pion correlations should be binned 
as differentially as possible. One may exploit longi¬ 
tudinal boost invariance and bin in the LCMS in 4D 
(A:T,gout,gside,®ong) as done in Ref. [2l|. 

In the case of partial coherence one has three quanti¬ 
ties which cannot be extracted from 2-pion correlations 
alone: One can, however, consider two ex¬ 

treme scenarios of coherent emission. In one case it is as¬ 
sumed that the space-time structure of coherent emission 
is identical to chaotic emission (ty = Ty). In another 
case one assumes that coherent emission is point-like at 
the center of the collision which might be expected for 
Bose-Einstein condensate (ty = 1). One may also con¬ 
sider additional scenarios with a suitable parametrization 
of tij (e.g. Gaussian structure with an assumed radius). 
The chaotic pair exchange magnitude, Ty, can then be 
extracted from Eq.[T3with various assumptions of G. Eor 
the purpose of studying partial coherence, built multi- 
pion correlation functions are defined with Eqs. [TSlfT^ 
but with Tij taken from lower order measurements. A 
minimization of the between built and measured cor¬ 
relations for each Qn bin can be used to estimate G for 
different assumptions of the coherent source profile. 


A. Extraction of Ty from C3 

In high-multiplicity collision events, femtoscopy lies in 
a clean region of phase-space where background correla¬ 
tions unrelated to QS and final-state-interactions (FSI) 
are negligible. However, in low-multiplicity events, back¬ 
ground correlations such as mini-jets are non- 

negligible and not exactly known. The extraction of Ty 
from 2-pion correlations can then be unreliable. Instead, 
one may exploit 3-pion cumulants as was done in Ref. 
for which 2-pion background correlations are explicitly 
removed. The 3-pion cumulant can be binned in 3D pair 
invariant relative momenta (< 712 , 923 , 931 )■ The pair ex¬ 
change magnitude can be parametrized with an Edge- 
worth or Laguerre expansion 

Ty =s£;w(i?9y)e-^'7o/2 (20) 

00 

E'w(R, 9y) = 1 + — Hn{Rqij) , (21) 

n=3 **qV 2j 

Ty =sLg(R9y)e-^7-*/2 (22) 

00 ^ 

Lg{Rqij) = 14- —^Ln{Rqij ). (23) 

n—1 

For the Edgeworth expansion, (Rqij ) characterizes de¬ 
viations from Gaussian behavior, Hn are the Hermite 
polynomials, and are the Edgeworth coefficients. For 
the Laguerre expansion, Lg{Rqij) characterizes devia¬ 
tions from exponential behavior, L„ are the Laguerre 
polynomials, and are the Laguerre coefficients. In both 
expansions, s is an additional scale parameter and R is 
the 2"'^ cumulant of the correlation function. For the case 
when Tw = I, i? is the standard deviation of a Gaussian 
source profile. For the case when Lg = 1, R is the FWHM 
of a Cauchy (Lorentzian) source profile. 

The 3-pion cumulant correlation can then be fit ac¬ 
cording to Eq. [T3] to extract the coefficients of the Ty 
parametrization. The resulting Ty can be used to build 
correlation functions. However, since 03 ( 912 , 923 , 931 ) is 
only a 3D projection of a 9D correlation function, there 
are somewhat more limitations to the accuracy of Ty ex¬ 
tracted from 3-pion cumulants than 2-pion correlations. 
An estimate for the bias on built correlation functions 
caused by limited dimensionality and binning will be dis¬ 
cussed in section X. 


V. RADII MEASUREMENTS 

Most often in experimental analyses, 2-pion QS corre¬ 
lations are used to extract information on the freeze-out 
space-time structure of the particle emitting source (e.g. 
radii). Less explored are the multi-pion QS correlations 
which may also be used to extract freeze-out radii (^.[2^ 
[^ . For such measurements, one typically assumes fully 
chaotic emission and parametrizes Ty in Eos. IlSlfTOl with 
the usual Gaussian or exponential forms. 
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Concerning the dimensionality of multi-pion projec¬ 
tions, often the correlations are projected and fit with 
ID Lorentz invariant relative momentum, Qn- However, 
a ID fit to the full correlation function is incorrect due 
to the summation of terms with different powers of Tij in 
Eqs.[T8]and[T9l Instead, one may isolate and fit only the 
3-pion cumulant, C 3 , for which the 2-pion symmetrization 
terms are explicitly removed. The cumulant may be fit 
in ID with a Gaussian function only since the definition 
of Q 3 coincides with the Gaussian relative momentum 
dependence of T 12 T 23 T 31 . For a spherically symmetric 
Gaussian source of radius R, Tij = e~ and there- 

fore T 12 T 23 T 31 = Q 3 / 2 , Por all other fit functions, it 

is more appropriate to fit in 3D ( 03 ( 512 , < 723 , <? 3 i))- In the 
case of 4-pion correlations, even a Gaussian fit to the cu¬ 
mulant correlation function is not correct since Q 4 does 
not coincide with the relative momentum dependence of 
TijTjkTkiTii. The source radius may be measured from 
C 4 through the 6 D pair q information by performing a fit 
according to Eq. [15] with the appropriate parametrization 
of the T factors. To avoid such high dimensional his¬ 
tograms, one may also compute Eq. [15] for each quadru¬ 
plet (online) for a variety of R choices and project against 
Q 4 . One may then perform a minimization offline to 
determine the best fit R value. 


A. Effect of partial coherence on extracted radii 


The extraction of femtoscopic radii in high-energy colli¬ 
sion data is universally done with a fully chaotic assump¬ 
tion of particle emission. It is well known, however, that 
the presence of partial coherence does not only suppress 
the intercept of QS correlations but may also modify its 
width Sii- With partial coherence in Eg. 1111 one may 
expect the same extracted radius only when Tij = tij 
and ^ij = 4 >ij (he. identical space-time structure of both 
chaotic and coherent components). Furthermore, the co¬ 
herent fraction must be momentum independent. In the 
more probable cases where these conditions are not satis¬ 
fied, partial coherence will modify the usual mapping of 
the correlation function width to the freeze-out radius. 
To illustrate this effect, 2-pion correlation functions are 
compared without coherence to two cases with partial 
coherence in Fig. |4| One observes that even with a sub¬ 
stantial coherent fraction of 20 %, the suppression at the 
intercept (g = 0) is very small. Moreover, above g = 0.15 
MeV/c there is actually an enhancement as compared to 
the fully chaotic case due to the quantum interference be¬ 
tween the chaotic and the smaller coherent source (wider 
correlation). Thus, in the presence of a smaller coherent 
component, the traditional fitting procedure will under¬ 
estimate the chaotic source radius. A single Gaussian fit 
to the red dash-dotted curve in Fig.|4|yields a radius 25% 
smaller than the chaotic source radius. 



FIG. 4. (Color online) Modification of 2-pion QS correlation 
functions in the presence of partial coherence. The correlation 
function for a fully chaotic spherically symmetric Gaussian 
source with R^h = 8 fm is shown with a solid black line. 
Dashed blue and dash-dotted red lines include an additional 
spherically symmetric coherent component with 7?coh = 1 fm 
for G — 0.2 and G = 0.5, respectively. 


VI. TREATMENT OF 
FINAL-STATE-INTERACTIONS 


The measurement of pure QS correlations is compli¬ 
cated by the presence of ESI (e.g. Goulomb repulsion). 
Fortunately, the Goulomb-l-strong wave-functions for 2- 
pion correlations are well known to high accuracy [^ . Isoll . 


X 




(24) 


where k* and r* are the momentum and relative sepa¬ 
ration evaluated in the PRF, 5c is the Goulomb s-wave 
phase shift, 77 = (fc*a)“^, ^ = k*r* -f fc*r*, p = k*r*, 
Ac = 27rr/le^'^'^ — 1 ]“^ is the Gamov factor, F is the con¬ 
fluent hypergeometric function, fc is the strong scatter¬ 
ing amplitude re-normalized by the long-range Goulomb 
forces, G is an s-wave Goulomb function, and a is the 
Bohr radius taking into account the sign of the interac¬ 
tion. The 2-pion FSI correlation can be computed by 
averaging the modulus square of the symmetrized wave- 
function over an assumed freeze-out source profile and 
dividing by the same average done with pure plane-waves 


(|V^-k»(r*) -h7Ak»(r*)P) 

(|e-*k*r* +e*k*r*|2) 


(25) 


Note that in the PRF, g = 2k*. Often, the averaging is 
done with a spherical Gaussian source profile. A more ac- 
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curate averaging would be one which includes the larger 
r* tails created by pions originating from resonance de¬ 
cays. However, due to the relatively large value of the 
pion Bohr radius (±388 fm for like/unlike sign pairs) as 
compared to the typical value of r* ~ 20 fm in high- 
energy heavy-ion collisions, the short-range structure of 
the source only mildly alters the functional form of K 2 {q). 
The averaging of K 2 {q) is usually confined to the core of 
particle production where r* < 100 fm. 


A. Multi-body FSI 


The exact wave-function of the n-body (n > 2) 
Coulomb scattering is unknown. However, solutions do 
exist in all asymptotic regions of phase-space . In par¬ 
ticular, the region of phase-space where all inter-particle 
spacings are sufficiently large is known as Hq. In the Hq 
asymptotic limit, the n-body Coulomb system is a sum 
of („^ 2 )! 2 ! non-interacting two-body systems [ 13 , HI] . It 
has been shown to be highly successful in describing ion¬ 
ization by electron impact (e, 2 e) in atomic physics [s^j . 
The relevant triplet kinetic energy for the applicability 
of the Hq ansatz was estimated in Ref. [s^l 


Etot > 0.2 


he 

i?(fm) ’ 


(26) 


where R is the estimated source radius. The characteris¬ 
tic Gaussian ID radii measured in pp/p-Pb and central 
Pb-Pb colliaipns at the LHC are about 1.5 and 9 fm, re¬ 
spectively [ 2 ^. The resulting threshold triplet energies 
are therefore 26 and 4 MeV for pp/p-Pb and Pb-Pb, 
respectively. For a symmetric triangular configuration 
of momenta in the triplet rest frame, the corresponding 
threshold for Q 3 is identical to the Etot value. Above 
the threshold, particles rapidly separate from each other 
such that the shorter range 3-body Coulomb forces (not 
treated in the Hq ansatz) become negligible. At the 4- 
pion level, the Hq ansatz should remain a good approxi¬ 
mation provided that each triplet Q 3 in the quadruplet is 
above the threshold. The 3-pion sym metrized Coulomb 
wave-function in Hq is given as in Ref. |3ll - l33j by 

V'pi,P2.P3(*1,^2,2:3) = 

_L[ei(9l2-12+-ji3-13-H923-23)/3,j,^2(,j,^2)<I>13(a:i3)4-23(a:23) 

V6 

+e49'l2^21+9l3^23-H9^3^13)/3<I,^2(a;2l)'I>13(a:23)4'23(2:i3) 

+ei(9'i2^32+9l3^3i-H9^3^2i)/3<j,^2(a;32)<I>l3(a:3l)4-23(a:2l) 

+ e49'l2-13+9;3-12-^9^3-32)/3<J,^2(3;i3)<I>13(a:i2)4-23(3;32) 

+ e49l2^23+9i3^21-r923^3l)/3<J,^2(a;23)'t>13(a:2l)4-23(3;3l) 
+e49'i2-31-H9l3-32+9A-12)/3,J,^2(x3l)<I-13(2:32)<I>23(*12)],(27) 


where qh = - pj, Xij = Xi - Xj, and ^ij{xki) = 

y/Ao{T]ij)F{-ir]ij, l,'i(k^jr*, ± k*^rl{)) is the Coulomb 
modulation factor and is evaluated in the ij PRF. The 
3-pion Coulomb correlation is then defined similar to K 2 
as 


^ 3 ( 912 , 913, 923 ) 


(l'^Pl,P2.P3 {XI,X2,X3)\'^) 
(l'0pi.P2.P3(a^i:a;2,a;3;4> = 1)|2) 


(28) 


In practice, the calculation of n-body Hq wave- 
functions is quite involved. Moreover, its application to 
experimental data requires a multi-dimensional calcula¬ 
tion in terms of the pair invariant relative momenta. For 
Coulomb correlations of 4*^ order, a 6 D calculation is 
required. As such a calculation is not very practical, a 
simplified approach is desirable. 

Another approach is the so-called generalized Riverside 
(GRS) method which treats the n-body FSI correlation 
as a product of 2-body FSI factors 

A3(qi2, gi3, 923) = K2{qi2)K2{qi3)K2(q2-i), (29) 

A4 (912,913, 914, 923, 924, 934 ) 

= K2(qi2)K2{qi3)K2{qiA)K2{q23)K2{q2A)K2{q3A). (30) 

The GRS approach easily allows a multi-dimensional es¬ 
timate of multi-body FSI since only pair calculations are 
required. Galculations of both approaches are presented 
in Sec. IIXI 


VII. TREATMENT OF DILUTION FROM 
NON-FEMTOSCOPIC SEPARATIONS 

The isolation of pure QS correlations is further compli¬ 
cated by the contribution of pions from long-lived emit¬ 
ters. Such pions are typically separated from other pions 
by many tens to hundreds of femtometers at freeze-out for 
which the expected quantum interference peak is only at 
extremely low 9 (< 5 MeV/c) [s^. In the moderate 9 re¬ 
gion, pions from long-lived emitters effectively dilute the 
correlation function. Long-lived emitters include weak- 
decays (secondary contamination) as well as long-lived 
resonance decays such as the ij. 

The treatment of the dilution for 2-pion correlations 
is usually done with the following formula connecting 
the measured and QS distributions/correlation functions 

M- 

N2{pi,P2) = (1 - /c)M(pi)Vi(P2) -I-/cF'2(9)V;^®(pi,P2)(31) 
C2(9) = (1 - fc) + flK2{q)C§^{q). (32) 

The distribution, N 2 {pi,P 2 ), can be experimentally 
formed by sampling both pions from the same event 
whereas Ni(pi)Ni(j> 2 ) is sampled from separate events. 
The correlated fraction of pairs for which an observable 
femtoscopic correlation (QS±FSI) exists is denote d by 
/^. In the “core/halo” picture of particle production [l7| . 
particles may originate from one of two sources with very 
different characteristic radii . Those originating from 
short-lived emitters are from the core of particle produc¬ 
tion and contain an observable QS correlation with other 
core particles. Those from long-lived emitters are from 
the halo and do not observably interact with any other 
particle. In such a picture, the fraction of particles from 
the core is denoted by fc- 

The treatment of dilutions in the core/halo picture can 
be easily extended to 3-pion correlations. For the same- 
event 3-pion distribution, V 3 (pi,p 2 ,P 3 ), there are four 
possible configurations. In one case, all three pions origi¬ 
nated from the core for which the probability is ■ In the 
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second case, two pions originated from the core for which 
the probability is 3^(1 — /c). In the third case, only 
one pion originated from the core with a probability of 
3(1 — /c)^/c- In the final case, zero pions originated from 
the core with a probability of (1 — /c)^- The measured 
3-pion same-event distribution can then be deco mpo sed 
into its constituent parts while also utilizing Eq. I All fo r 
the decomposition of iV 2 (pi,P 2 ) into its two parts | 21 |. 
The same-event 3-pion distribution can be written as 

N3{pi,p2,P3) = f3lNl{pi)Nl{p2)Nl{p3) 

+ f32lN2(pi,P2}Ni{p3) + N2(p2,P3)Nl{pi) 

+ A^2(P3,Pi)A^i(P2)] 

+ [pi,P2,P3), ( 33 ) 

where is the quantity of interest to be extracted 
from the measured distributions (no QS superscript) and 
A'a is the 3-body ESI factor. The additional distribu¬ 
tion, N 2 {pi,Pj)Ni(jik), provides information from 2 -pion 
symmetrizations alone in the dilute gas limit. It may 
be experimentally measured by taking two pions from a 
single-event while taking the third from a mixed-event. 
In the core/halo picture, /31 = (1 - fc)^ + 3/^(1 - fc)^ - 
3(1 - U){1 - /2), /32 = 1 - /c, and /33 = 

Eour-pion correlations have four sequences of sym¬ 
metrizations as illustrated in Eig. [H Isolation of each 
sequence can be accomplished with the following addi¬ 
tional distributions: 


Ni{p3)N3{p2)Nr{p3)Nr{pi), (34) 

N2ipi,P2)Ni{p3)Ni{p4), (35) 

N2{pi,P2)N2{P3,P4), (36) 

N3{pi,P2,P3)Ni{p4), (37) 

^4(P1,P2,P3,P4), (38) 


which can be experimentally measured by taking the ap¬ 
propriate number of particles from t he s ame event and 
the rest from mixed events. Eor Eq. [36l one takes two 
particles from one event and the other two from a sepa¬ 
rate event. Extending the treatment of dilutions to the 
4-pion case, one obtains 

N4{pi,P2,P3,P4) = f4lNl(pi)Ni(p2)Ni{p3)Ni{p4) 

+ f42W2(pi,P2)Ni{p3)Ni{p4) 

+ ‘13’ -I- ‘14’ -I- ‘23’ -I- ‘24’ -I- ‘34’] 

-t f43[N3{pi,P2,P3)N4{p4) 

+ ‘124’ -I- ‘134’ -I- ‘234’] 

+ f44K4N^^{pi,P2,P3,P4), (39) 

where is again the quantity of interest to be ex¬ 

tracted from the measured distributions and the 4-body 
ESI correlation, K 4 . Permutations indicated by Hj’ rep¬ 
resent N 2 {pi,Pj)Ni{pk)Ni{pi). Permutations indicated 
by Hjk' represent N 3 {pi,pj,pk)Ni{pi). In the core/halo 
picture, /41 = -3(l-/c)'‘-8/c(I-/c)^+6(I-/c)(l-/c)^ 
/42 = -(I - fc)^, /43 = (1 - /c), /44 = Another pos¬ 
sibility not represented in Eq. [39] is given by two pairs of 
interacting pions, N 2 {pi,P 2 )N 2 ip 3 ,P 4 )- This possibility 
is absent in the core/halo picture and will be shown to 
be quite small in a more realistic model in Sec. El 


A. A parameter 

Most often in experimental analyses, the parameter 

in Eq. [32] is replaced with the so-called A param- 
eter and is parametrized without coherence as: 

C 2 ^ = 1 -|- [36l - [3^ . Such a parametrization 

can be flawed in regards to ESI corrections (AA' 2 ) in two 
noteworthy ways. Eirst, when the true functional form 
of C 2 is not known, a mismatch between the true and 
assumed function form can cause a large bias on the ex¬ 
tracted A parameter M- As the ESI correction is the 
product of A and A 2 , a biased value of A also biases the 
ESI correction. Second, coherent emission leads to the 
suppression of but not of ESI correlations. In both 
cases, a single suppression parameter (A) is not appropri¬ 
ate. The problem can be circumvented by parametrizing 
the QS correlations instead by = 1-1- , 

where in addition to in Eq. |32l a second suppression 
parameter is used. 

As in Ref. [2l| , may be determined less ambiguously 
from a fit to tt+tt” correlation functions for which only 
ESI correlations contribute. Assuming the equivalence 
of /+ with f~ as expected at high energies Q, tt+tt" 
correlations can be fit according to C 2 = (1 ~ /c) + 
fcK 2 {q,R)- Owing to the large value of the pion Bohr 
radius (388 fm) as compared to the typical relative sep¬ 
aration at freeze-out in high-energy heavy-ion collisions 
(~ 20 fm), 7r+7r“ ESI correlations are much less sensi¬ 
tive to the short-range structure of the source. A fit to 
7r+7r“ correlations is then simpler in practice as K 2 {q, R) 
changes less rapidly for different source profiles with the 
same characteristic radius. 

The chaotic upper limit to the conventional A param¬ 
eter is unity. However, experiments overwhelmingly re¬ 
port values less than unity-sometimes as low as 0.3 [ 3 ^ . 
There can be several sources of suppression for the con¬ 
ventional A parameter which are listed below in approx¬ 
imate order of magnitude: 

• Dilution from pions of long-lived emitters 

• Gaussian fits to non-Gaussian correlations 

• Detector momentum smearing 

• Pion mis-identification 

• Einite q bin width 

An additional more interesting source of suppression not 
listed above, is the possibility of pion coherence. How¬ 
ever, the suppression is expected to be rather small for 
2-pion correlations which was estimated to be ~ 0.05 at 
the LHG [m . The exact order of the above contributions 
depends on the experimental conditions of the analysis. 
One of the main reasons for low values of A is often due 
to the use of Gaussian fits to correlation functions which 
are intrinsically non-Gaussian. In Ref. [2l|, an estimate 
can be made for the top three sources of suppression. A 
Gaussian fit to C'^''’(g) yielded A ^ 0.4. However, a fit 


to C 2 {q) gave A ~ 0.7, which indicates a suppression of 
about 0.3 from long-lived emitters in ALICE. The sepa¬ 
rate analysis of revealed a possible coherent fraction 
of about 23% for which the suppression to the intercept 
of 2-pion correlations is 0.23^ = 0.05. This leaves a re¬ 
maining suppression of about 0.25 due to non-Gaussian 
features in the correlation function. 

Note that the non-Gaussian features for a ID analysis 
are expected to differ from 3D analyses (gout, 9side, ®ong) 
even if the source is Gaussian in all three dimensions 
but with different radii in each dimension. The STAR 
collaboration reported values between 0.3 and 0.4 for 3D 
Gaussian fits which likely indicates substantial non- 
Gaussian features even in 3D heavy-ion analyses. 


VIII. CUMULANT AND PARTIAL CUMULANT 
CORRELATION FUNCTIONS 

The full n-pion correlation function contains the full 
set of symmetrizations as previously illustrated. One 
may isolate different levels of symmetrization by sub¬ 
tracting different types of n-pion spectra. The 3-pion 
cumulant correlation function can be defined as 

C 3 = [N^^{P1,P2:P3) - 6 iN^^{pi,P 2 )Ni{p 3 ) 

+ 62 Nf]/Nl (40) 

where Nf denotes Ni{pi)Ni{p 2 )Ni{p 3 ) for brevity. For 
the case of same-charge triplets, 61 = 62 = 3. For the 
mixed-charge case, i5i = i 52 = 1. The cumulant correla¬ 
tion function removes the 2 -pion symmetrizations. 

With 4-pion correlations, two types of partial cumu- 
lants as well as the full cumulant correlation function are 
defined: 

a4 = [n2^{pi,P 2,P3,P4) - eiN^^{pi,p2)Nl{p3)Nl{p4) 

+ eiNf]/Nf, ( 41 ) 

b4 = W2^{Pl,P2,P3,P4) -3Af^®(pi,P2)N^®(P3,P4) 

+ 3Nf]/Nf, ( 42 ) 

C4 = [N^®(pi,P2,P3,P4) - C2N^®(Pl,P2,P3)Nl(p4) 

- <^3^^^ {pi,p2)N^^ {p3,p4,) + {pi,p2)Nl{p3)Nl{p4,) 

+ {e2+e3-e4)Nf]/Nf. ( 43 ) 

The cumulant and partial cumulants isolate different 

sequences of symmetrization. When a mixture of iden¬ 
tical and non-identical pions is used (e.g. Tr"*" and 7 r“ 
mixtures), one need only subtract the symmetrizations 
from identical pion pairs and triplets. For same-charge 
quadruplets the coefficients are: ei = 6 ,e 2 = 4, eg = 
3 ,e 4 = 12. For mixed-charge quadruplets of type 1 
(T ± ±±) the coefficients are: ei = 3, 62 = 1, 63 = 0, £4 = 
0. For mixed-charge quadruplets of type 2 (=F T ii) the 
coefficients are: ei = 2,62 = 0, £3 = I, £4 = 0. 

In the case of same-charge quadruplets, the first partial 
cumulant, a 4 , is defined such that all six terms of 2 -pion 
symmetrization are removed. The second partial cumu¬ 
lant, b 4 , further removes the three 2 -pair symmetrization 


terms. Note that contains two terms of 2-pion sym¬ 
metrization as well as one term of 2 -pair symmetrization. 
Finally, the cumulant further removes the eight terms 
of 3-pion symmetrization and is thus a measure of the 
genuine 4-pion symmetrization alone, Note that V 3 con¬ 
tains six terms of 2 -pion symmetrization and two terms 
of 3-pion symmetrization. The above equations for the 
cumulant and partial cumulants are valid for identical 
pions. When a mixture of identical and non-identical pi¬ 
ons is used (e.g. 7 r+ and 7 r“ mixtures), one need only 
subtract the symmetrizations from identical pion pairs 
and triplets. 

The cumulant and partial cumulant correlation func¬ 
tions have an advantage over the full correlation functions 
(Eq. [T]) in low multiplicity collision events as pointed out 
in Ref. [ 2 ^. In such events, the contributions from non- 
Bose-Einstein correlations (e.g. mini-jets) have a non- 
negligible effect. Two-pion background correlations are 
explicitly removed with the cumulant and partial cumu¬ 
lants. The remaining higher order background correla¬ 
tions were shown to be very small in Ref. [IJ . Also note 
that the choice of the Q 3 projection variable may also 
influence the flatness of the baseline. QS correlations are 
localized at low Qa while 3-pion mini-jet correlations may 
be less well localized with such a variable (smeared). 


1. Practicality: Nested loops 


In high multiplicity events such as those produced in 
heavy-ion collisions, the analysis of multi-particle correla¬ 
tions directly with nested loops is prohibitively expensive 
in terms of GPU time. One can circumvent this problem 
by only retaining pion-pairs in a specific region of q. A 
system of switches using 2D arrays of Boolean variables 
can be formed for each type of pion pair (e.g. pion-I 
from event A and pion-2 from event B). Two classes of 
switch-arrays should be formed: one for the femtoscopic 
analysis region at low q and one for the normalization 
region at high q. The cutoff for the low q region should 
be chosen to include the dominant region of QS corre¬ 
lations. The width of the normalization region should 
be chosen to contain sufficient statistics while still mini¬ 
mizing CPU time. Nested loops are still utilized for the 
full analysis except that at the start of 2 "'^ and higher 
loops one checks each pair with the switch-array. If 
any pair is turned “off”, one skips the entire n-tuple. 
In this manner, undesired n-tuples are skipped before 
any time consuming calculations are performed. The 
algorithm is depicted with the following pseudo-code. 
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Stage one: Set pair Boolean switch 
for all z = 0 to TV do 
for all j = z + 1 to do 

If Qmin ^ Q ^ 9max then 

SWITCH H[j] = 1 
else 

S WITCH = 0 

Stage two: Perform 4-pion correlation analysis 

for all z = 0 to TV do 

for all j = z + 1 to do 

if SWITCH[i][j] == 0 then 

skip 

for all /c = j' + 1 to A^ do 
if SWITCH[i][fc] == 0 then 
skip 

if SWITCH[j][fc] == 0 then 
skip 

for all ^ = fc + 1 to A^ do 
if SWITCHH[;] == 0 then 
skip 

if SWITCH[j][/] == 0 then 
skip 

if SWITCH[A:][Z] == 0 then 
skip 

Store quadruplet 
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FIG. 5. (Color online) Triplet fractions in therminator. 
Blue triangles (3 interacting pairs), red squares (1 pair), green 
triangles (2 pairs), and black circles (0 pairs) are shown. The 
core/halo values (given = 0.803) are shown with blue (3 
pairs), red (1 pair), and black (0 pairs) lines. 


IX. THERMINATOR CALCULATIONS 

In this section several calculations of 3- and 4-pion 
correlations in the heavy-ion generator therminator 2 
[3^ [4^ are provided, therminator 2 is a Monte Carlo 
event generator based on the statistical hadronization of 
partons at kinetic freeze-out. A small number of physi¬ 
cal input parameters are used such as the temperature, 
chemical potentials, initial size, and the velocity of col¬ 
lective flow which has been shown to describe a great 
deal of the RHIC and LHC observables. An important 
feature of therminator is the inclusion of the full set 
of hadronic resonances. 

A 3-l-lD viscous hydrodynamic hypersurface is input 
into the statistical hadronization process. Events were 
generated with the following settings: Pb-Pb collisions 
at y^SNN = 2.76 TeV, b = 2.3 fm (impact parameter), 
Ti = 512 MeV (initial central temperature), to = 0.6 fm 
(starting time of hydrodynamics), 7/ = 140 MeV (freeze- 
out temperature). Approximately 270 x 10^ events were 
generated. Only charged pions were retained with the 
following criteria: 0.16 < pt < 1-0 GeV/c, Izy] < 0.8. To 
simulate the typical experimental resolution of a track’s 
distance-of-closest approach (DCA) to the primary ver¬ 
tex, pions which freeze-out greater than 1 cm from the 
collision center are rejected. For the correlation part of 
this analysis, only the dominant region of QS correla¬ 
tions are analyzed by requiring a cut on the pair relative 
separation in the PRF: 0.1 < r* < 100 fm. 


A. / coefficients 

The / coefficients in Eqs. [33] and |39| characterize the 
probabilities of short- and long-range interactions for 
triplets and quadruplets. A pair is deemed interacting 
in THERMINATOR if its value of r* is less than a cer¬ 
tain cutoff, which was taken to be 80 fm (well above the 
dominant QS region). With this cutoff, the fraction of 
interacting pairs is = 0.803. The triplet fractions are 
shown in Fig. [5] The blue, red, and black lines correspond 
to the core/halo values of /as, 3 / 32 , and / 31 , respectively. 
Note that the case of 2 interacting pairs is absent in the 
core/halo picture but non-zero in general. One example 
of such a case would be a linear configuration of 3 pi¬ 
ons where the middle pion is separated from both end 
pions by an amount less than the cutoff while the end pi¬ 
ons themselves are separated by a greater distance. The 
quadruplet fractions are shown in Fig.| 6 | The black, blue, 
red, and green lines correspond to the core/halo values 
of / 44 , 4 / 43 , 6 / 42 , and /41 respectively. Note that the 
cases of 5, 4, and 2 interacting pairs are absent in the 
core/halo picture but non-zero in general. The case of 
5 interacting pairs can occur for a parallelogram config¬ 
uration of the 4 pions where the longest axis represents 
the length above the r* cut. An example of a case with 4 
interacting pairs can occur with an equilateral triangle of 
3 pions and the fourth pion being located outside of the 
triangle but near one of the vertices. The case of 2 in¬ 
teracting pions can a be a linear configuration of 3 pions 
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FIG. 6 . (Color online) Quadruplet fractions in therminator. 
From top to bottom, the markers indicate 6 , 3, 1, 5, 4, 2, and 
0 interacting pair(s). The core/halo values (given = 0.803) 
are shown with black (6 interacting pairs), blue (3 pairs), red 
(1 pair), and green (0 pairs) lines. 
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Triplets 

hi 

fi2 

hi 

r* < 60 

0.675(4-9%) 

0.094(-17%) 

0.041(-25%) 

r* < 80 

0.720(4-6%) 

0.083(-15%) 

0.030(-24%) 

r* < 100 

0.738(4-6%) 

0.079(-15%) 

0.026(-24%) 


TABLE I. Triplet / factors in the core/halo picture and the 
THERMINATOR modification percentage in parenthesis. Note 
that 3/32 is shown in Fig. [5l 


Quadruplets 

/44 

hi 

/42 

/41 

r* < 60 

0.593(4-10%) 

0.083(-12%) 

0.012(-26%) 

0.007(-18%) 

r* < 80 

0.645(4-6%) 

0.075(-9%) 

0.009(-19%) 

0.004(-ll%) 

r* < 100 

0.667^5%) 

0.071(-8%) 

0.008(-18%) 

0.003(-7%) 


TABLE 11. Quadruplet / factors in the core/halo picture 
and the therminator modification percentage in parenthe¬ 
sis. Note that 4/43 and 6/42 are shown in Fig. [O] 
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FIG. 7. (Golor online) C 2 versus q calculated in thermina¬ 
tor. A Gaussian as well as Edgeworth fit is shown. The fit 
range is that which is shown in the figure. 0.2 < Ict < 0.3 
GeV/c. 


and the fourth being far away from the rest. For both 
triplet and quadruplet fractions there is no significant Qa 
and <54 dependence. Thus, the / coefficients are largely 
independent of relative momentum. 

In experiment, one cannot isolate each of the possibil¬ 
ities with mixed-event techniques. For instance, the case 
of 2 interacting pairs in the triplet cannot be isolated 
using the three available distributions: N^{pi,p 2 iP 3 ), 
N 2 {pi,P 2 )Ni{pz), and Ni{pi)Ni{p 2 )Ni{pz). The mod¬ 
ification to the core/halo / coefficients in Eqs. [33] and 
ISniare estimated by combining the appropriate thermi¬ 
nator fraction types. For triplets, the modification to 
/as is given by the sum of 2 and 3 interacting pair frac¬ 
tions in Fig. |5l For quadruplets, the modification to /44 
is given by the sum of 4, 5, and 6 interacting pair frac¬ 
tions in Fig. ini The modification to /43 is given by the 
sum of 2 and 3 interacting pair fractions. Table |I| shows 
the core/halo values of /as, fz 2 ^ /ai and their percentage 
modification in therminator in parenthesis. The calcu¬ 
lation is done for three different choices of the r* cutoff. 
Table El shows the core/halo values of / 44 ,/la,/ 42,/41 


and their percentage modifications in therminator in 
parenthesis. 


B. Measured and built correlation functions 

The 2-, 3-, and 4-pion QS correlation functions in 
THERMINATOR are now presented. The correlation func¬ 
tions are calculated via the fully symmetrized 2-, 3-, and 
4-pion plane-wave functions. The symmetrized 2-pion 
plane-wave function is given by'l '2 = ]■ 

The symmetrized 3-pion plane-wave function was given 
in Eq. |27l with $ = 1. In appendix A the sym¬ 
metrized 4-pion plane-wave function is given. The mod¬ 
ulus square of the plane-wave functions is applied as a 
pair/triplet/quadruplet weight and is averaged over the 
freeze-out space-time coordinates in each event. The 
weight is then averaged over all events. The 2-pion 
correlation function is shown in Eig. E It is fit with 
a Gaussian and Edgeworth parametrization: C 2 {q) = 
I -I- s'^E^{Rq) e~^ . The Gaussian case corresponds 
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FIG. 8. (Color online) Tliree-pion correlations versus Q 3 cal¬ 
culated in THERMINATOR. A Gaussian as well as Edgeworth 
fit to the cumulant correlation function is shown. The fits are 
performed in 3D ( 512 , 923 , 531 ). The built correlation func¬ 
tions are also shown as block histograms. 0.16 < At ,3 < 0.3 
GeV/c. 

to Aw = 1- Note, that only femtoscopically separated 
particles are retained (r* < 100 fm) and FSI were not 
included in the weighting procedure. In our calculation, 
due to the full symmetrization, pion production is fully 
chaotic and thus the s parameter of the fit is just a scale 
factor unrelated to coherence. 

The 3-pion correlation function projected against Q 3 
is shown in Fig. HI A Gaussian and Edgeworth fit to the 
cumulant correlation is performed in 3D according to 

C3(512,523,531) = 1-t ^'^12''"'^23''"'^31^^^Fw,i2i?w,23Fw,31, 

(44) 

which can be obtained by substituting Edgeworth ex¬ 
change amplitudes from Eq. HD] into Eq. |TD| (neglect¬ 
ing the 3-pion phase). Two types of built correlation 
functions are shown in Fig. HI Dashed blue lines cor¬ 
respond to C 3 built from extracted directly from 
C 2 (fcT, 9out, <Zside, ®ong)- Solid blue lines correspond to 
C 3 built from extracted from 03 ( 512 , 523 , 931 ) fits. 
The building of correlation functions is done according 
to Eq. [TH] with G = 0. The ratio of measured to built 
C 3 is shown in Fig. jD] Built correlations are found to re¬ 
produce the measured ones to within a few percent while 
those from G 2 (A:t, 5out, 5side, 5iong) are found to be some¬ 
what more accurate. This is partly expected since 2-pion 
correlations are binned much more differentially. 

The 4-pion correlation function projected against Qd 
is shown in Fig. 1101 The partial cumulants, a 4 and b 4 
represent the removal of 2 -pion and 2 -pion-|- 2 -pair sym- 
metrizations, respectively. The full cumulant, C4, rep- 


FIG. 9. (Color online) Ratio of measured to built 3-pion cor¬ 
relation functions. 



FIG. 10. (Golor online) Four-pion correlations versus Q 4 cal¬ 
culated in THERMINATOR. The built correlation functions are 
also shown as a block histogram. 0.16 < At ,4 < 0.3 GeV/c. 


resents the isolation of genuine 4-pion symmetrizations. 
Built correlation functions are shown with the colored 
box histograms. The building of correlation functions 
is done according to Eq. [TD| with G = 0. The ratio of 
measured to built G 4 is shown in Fig. [TTJ Similar to 
built 3-pion correlations, built 4-pion correlations gener- 
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(GeV/c) 



FIG. 11. (Color online) Ratio of measured to built 4-pion 
correlation functions. 


ally under predict the measured ones by a few percent. 
The reason for the slight bias is investigated further in 
Sec. IIX Cl with ra and is related to the finite binning 
of multi-dimensional correlation functions. For the built 
correlations in this subsection, 2-pion correlations are 
binned as follows: 4 /ct bins dividing the full range, 
5 MeV/c bin widths for (7out,9side,®ong- Linear interpo¬ 
lation is used in the dimension while cubic is used in 
the other dimensions. Three-pion cumulants are binned 
with 5 MeV/ c widths in each pair . 


1. Including partial coherence 

To demonstrate the building of QS correlations with 
partial coherence, a coherent source of pions is inserted 
in THERMINATOR. Pions are randomly reassigned such 
that 35% of the pions in an event originate from a spher¬ 
ical Gaussian coherent source with i?coh = 1 fiR- That 
is, G = 0.35 and tij = e“ The momentum dis¬ 

tribution of the coherent pions is identical to that of the 
chaotic pions. The expected suppression from coherence 
is incorporated by setting the appropriate plane-wave 
contributions to zero. In the fully symmetrized plane- 
wave function, individual plane-waves corresponding to 
an n-pion symmetrization are set to zero if more than 
n — 1 of the pions are from the coherent source. 

The 2-pion correlation function with the parametrized 
coherent source is shown in Fig. [121 The 3-pion corre¬ 
lation function with the parametrized coherent source is 
shown in Fig. 1131 An Edgeworth fit to the cumulant cor¬ 
relation is shown with a dashed blue line. The fit well 


FIG. 12. (Golor online) C 2 versus q calculated in thermina- 
TOR with G = 0.35 and Uj = 913 /2 ^ Gaussian as well 

as Edgeworth fit is shown. Further details are the same as in 

Fig.m 



63 (GeV/c) 


FIG. 13. (Golor online) Three-pion correlations versus Qs cal¬ 
culated in THERMINATOR with G = 0.35 and Uj = 

An Edgeworth fit is shown. Further details are the same as 
in Fig. m 
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FIG. 14. (Color online) Ratio of measured to built 3-pion 
correlation functions with G = 0.35 and tij = e“ 

describes the cumulant for Q 3 > 0.015 GeV/c. Below 
that value, the fit is below the data and is part of the 
reason why the built correlation from C 3 underestimates 
the measured C 3 in the same region. The built correla¬ 
tion function from C 2 is also beneath the measured C 3 
for reasons described in Sec. lIX Cl The ratio of measured 
to built C 3 is shown in Fig. [TH The 4-pion correlation 
function with the parametrized coherent source is shown 
in Fig. [TS] The ratio of measured to built is shown in 
Fig. [161 The built correlations incorporate the known G 
and tij values and thus illustrate how built correlations 
can be used to estimate the coherent fraction. 


C. rs and r4 measurements 

The 3- and 4-pion phase factors in therminator are 
measured with ra and in Eqs. [T3] and 1161 In the 
computation of both denominators, the pair exchange 
magnitudes appear. As this quantity must be evaluated 
numerically by averaging over the sources produced in 
many events, it is tabulated discretely in the first pass 
over the data. In the second pass one interpolates be¬ 
tween the bins linearly or cubically. Finite bin width is 
found to cause a bias on ra and r 4 in the same manner as 
it occurred for the built correlation functions. The bias 
is largest for linear interpolation between widely spaced 
bins of (?out j 9 side, and giong and depends on the concav¬ 
ity of the correlation function. For instance, negative 
concavity with linear interpolation generally leads to an 
underestimation. The concavity of a Gaussian correla¬ 
tion function is negative for q < y^lj 2 R'^ and positive 


FIG. 15. (Golor online) Four-pion correlations versus Q 4 cal¬ 
culated in THERMINATOR with G — 0.35 and tij = 

Further details are the same as in Fig. 1101 



FIG. 16. (Color online) Ratio of measured to built 4-pion 

d 2 2 / o 

correlation functions with G — 0.35 and tij = e“ ' . 

for larger q. The bias may be reduced somewhat with a 
cubic interpolator. 

The 2-pion Tjj factors are computed in 3 different ways 
depending on the number of q bins and on the interpola¬ 
tion scheme between q bins (linear or cubic). The de¬ 
pendence of Tij was found to be fairly linear and thus lin- 
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FIG. 17. (Color online) rs calculated in therminator with 
various types of Tij tabulations/interpolations. The quartic 
fit is applied to the solid points. Statistical errors are only 
drawn for the solid points and are the same for the others. 
The chaotic upper limit is shown with the red dotted line. 


FIG. 18. (Color online) r 4 calculated in therminator with 
various types of Tij tabulations/interpolations. The chaotic 
upper limit is shown with the red dotted line. 


ear interpolation is always used between /cx bins. Three 
different q bin widths (1, 2, or 5 MeV/c) are tried while 
dividing the full interval into 4 bins. In Fig. [IT] and 
[iHl rs and r 4 calculations in therminator are shown. 
Similar to Ref. [ 2 l|, quartic {r„{Qn) = — oQn)) 

are shown for both and r 4 for Aq = 5 MeV/c, cubic. 
The quartic intercept parameter to ra is / = 2.08 ± O.OI 
while a is consistent with zero. For r 4 they are / = 
6.6 ± 0 . 1,0 = (6 ± 2) X 10^. The near independence of 
ra and r 4 with Qa and Q 4 indicate a negligible effect of 
3- and 4-pion phases in therminator. This feature is 
expected for the case where the space-time point of max¬ 
imum pion emission is momentum independent fl8l |. In 
Fig. [13 one also observes that the chaotic upper limit 
(2.0) is exceeded to varying degrees. The bias is caused 
by the limited dimensionality and finite bin widths of 
Tij, especially in the case of linear interpolation. Cubic 
interpolation is observed to reduce the bias. 

With a coherent component inserted as in Sec. IIXB 1 [ 
the expected suppression of ra is observed as seen in 
Fig. [ini The quartic fit parameters in Fig. [12] are: 
I = 1.77 ± 0.01 and a = (1.6 ± 0.5) x 10^. The inter¬ 
cept parameter is related to the coherent fraction as 


I = 2VT^ 


1 + 2G 
(l + G)3/2’ 


(45) 


which yields G = 0.33±0.01. Note that in the recent AL¬ 
ICE ra measurement [2lj, linear interpolation was used 
with 4 fcx bins and 5 MeV/c q bin widths. Correction 
for the finite binning effect may lower ra and increase 


the extracted coherent fraction. Concerning the a pa¬ 
rameter, it was shown in Ref. that quartic behavior 
is expected when the space-time point of maximum pion 
emission is momentum dependent. As the a parameter 
for the chaotic part of the therminator source is zero 
in Fig. 1171 and since the momentum spectrum of the co¬ 
herent and chaotic components are identical, one expects 
a to vanish also with the inserted coherent component. 
The non-vanishing value of a also characterizes the same 
bias as attributed to /. 


D. 3 -pion FSI calculations 


In Ref. [2l| and Fig. [201 a comparison between Hq and 
GRS 3-pion FSI correlation is made. One observes a 
striking similarity between the two methods, especially 
for same-charge triplets. Less similar is the mixed-charge 
calculation for which the cumulants differ by as much as 
30%. In Ref. [2l|, the mixed-charge calculations were 
more similar and is due to the inclusion of strong FSI 
in the factorization ansatz. The calculation in Fig. [201 
excludes strong FSI. The similarity of the GRS and Oq 
method is related to the factorizability of separate F fac¬ 
tors as well as that between F and plane-wave factors. 
As one decreases the source size in the calculation, both 
methods converge as the F factors approach unity. 
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FIG. 19. (Color online) rs calculated in therminator with 
G — 0.35 and tij = " 21^/2 gjjj widths are 5 MeV/c with 

cubic interpolation. A quartic fit is shown. Statistical errors 
are only drawn for the solid points and are the same for the 
others. 


X. MULTI-BOSON DISTORTIONS 


The standard framework of Bose-Einstein interferom¬ 
etry neglects the effect of higher-order symmetrizations 
and is only valid in the limit of low phase-space densities 
at freeze-out [iMl- In such a dilute pion gas, n-pion 
correlation functions are calculable purely in terms of 
n-pion symmetrizations. At higher densities, the higher- 
order symmetrizations become significant and distort the 
correlation functions. Previous calculations revealed a 
widened 2-pion correlation function as well as a sup¬ 
pressed intercept, both of which are unrelated to quan¬ 
tum coherence. Here, the effect of the distortions on the 
comparison of built and measured 3- and 4-pion correla¬ 
tion functions is analyzed. 

The computation of all higher orders of symmetriza¬ 
tions is greatly simplified when the n-pion emission func¬ 
tion is assumed to factorize into a product of single parti¬ 
cle emission functions S{x^p). A remaining complication 
is the computation of the n\ plane-wave functions which 
caj b& alleviated with the help of the Pratt ring integrals 



FIG. 20. (Color online) Comparison of same and mixed- 
charge three-pion FSI correlations. Hq and generalized River¬ 
side (GRS) methods are shown. The calculation was per¬ 
formed in THERMINATOR. The bottom panel shows the differ¬ 
ence between the two methods, AK 3 — K 3 {flo) — Ar 3 (GRS), 
divided by As (GRS) - 1. 


spectrum at a fixed multiplicity n are given in Ref. 


n —1 

M^’(P) = E 


l{n - 1 - j)\ 

UJnln\ 


G, 


0 + 1 ) 


(p. p) 


N,i^^(pi,P2) 


_ ‘^n — 2—j / {ri — 2 — j)\ 

3 

X X^Gi+l(pi,pi)Gj_i+l(p2,P2) 

i=0 

+ Gi+i(pi,p2)Gj_i+i(p2,pi) 


(48) 


(49) 


where Wn is the Bose-Einstein weight of an event with n 
identical bosons. Extending the techniques in Ref. [l^ . 
the order i spectra are obtained by the appropriate sym- 
metrization oi i G factors. The 3-pion spectrum is given 
by 


Gi(pi,P2) = J d‘^xSi^(pi+p2),x)e^<-P^-P^'>^, (46) 

G„(pi,P 2) = y’d3k2...d3k„Gi(pi,k2)...Gi(k„,p2). (47) 

With a Gaussian ansatz for the single particle emis¬ 
sion hmction, the ring integrals are given analytically in 
Ref. [13 and are used here to estimate the distortion to 3- 
and 4-pion correlation functions. The single and 2-pion 


M®'(PRP2,P3) = 

2=0 


n — 3 


‘^71—3—jI{n — 3 — j)\ ^ 

1=0 m=0 


Ldnlnl 


'y ^ G;_)_i (pi, Pcti )G^m+l (p2 ) Pq'2 '}^j — l — 7n-\-l (P 3 ? Po'3 ); (^0) 

q:C3) 


where the set of all permutations of the 3 pions is given 
by and 01 . 2,3 represent the permuted indices. The 
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4-pioii spectrum is given by 


M'‘\pi,P2,P3,P4) = ^ 


u)n-i-jl{n - 4 - j)! 




ui„ln\ 


y ] Gz + l(pi 1 Pai iGrn-l-l (P2 , Pa2 ) 




^ G^s + 1 (p3 ! Pq ;3 )Gj — l — m — s-\-l (P4 ; Pq ;4 )? 


EE E 


(51) 


where the set of all permutations of the 4 pions is given 
by and ai, 2 , 3,4 represent the permuted indices. The 
normalization constant for the correlation functions are 
obtained by the term with the lowest order ring integrals 
(Gi). The 3- and 4-pion normalizations are 


(n 


{n - l)(n - 2 )ujluJn-3 


l)(n - 2)(n - 3 )uj^uJn-4 


(52) 

(53) 


In the Gaussian model of Ref. [l^, three parameters 
characterize the pion spectra: A, rg, and n. A charac¬ 
terizes the width of the momentum spectrum unaffected 
by the Bose-Einstein enhancement. The isotropic spa¬ 
tial width of the emission function in the PRF is given 
by tq. To quantify the multi-boson distortions in the re¬ 
cent LHC central 0 — 5% Pb-Pb data [2l|,[2j, A = 0.25 
GeV/c, ro = 8.5 fm, n = 700,1400, and (p) = 0.2 GeV/c. 
The choice of rg is based on the extracted radii from ID 
Gaussian fits to the measured correlation functions [^ . 

In Fig. [2T]|23l the distorted and undistorted 2-, 3-, 
and 4-pion correlation functions are compared. It is 
clear that the multi-boson distortions widen all correla¬ 
tion functions while there is no significant suppression of 
the intercept, which occurs for much higher pion densi¬ 
ties. An important feature is the similarity of built and 
distorted correlation functions. This is due to the fact 
that the 2-pion correlation, which is used in the building 
procedure, is distorted as well. The ratio of distorted 
to built 3- and 4-pion correlation functions stays within 
O.OI below unity. The calculations were also performed 
with rg = 7 fm and n = 700 which result in larger pion 
densities. The increased widening of the 2-pion corre¬ 
lation function leads to a 10% underestimation of rg, 
when extracted with a Gaussian fit function. The ratio 
of distorted to built 3- and 4-pion correlation functions 
remained very similar. The proximity of the ratio with 
unity demonstrates the robustness of the built approach 
to multi-pion distortions at moderate pion densities. 


XI. SUMMARY 

Techniques to isolate and analyze multi-pion QS corre¬ 
lations have been presented. The concept of constructed 
or built correlation functions may be used to estimate the 
coherent fraction at measured values of relative momen¬ 
tum instead of at the unmeasured intercept of correlation 
functions. 



q (GeV/c) 


FIG. 21. (Color online) Multi-boson distorted 2-pion correla¬ 
tion function compared to the undistorted case. The setting 
with n = 700 roughly correspond to 0 — 5% Pb-Pb collisions 
at the LHC. 



^3 (GeV/c) 


FIG. 22. (Color online) The multi-boson distorted 3-pion 
correlation function compared to the undistorted case. Also 
shown is the built correlation function constructed from the 
distorted 2-pion correlation function in Fig. [21] The setting 
with n = 700 roughly correspond to 0 — 5% Pb-Pb collisions 
at the LHC. 
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(GeV/c) 


FIG. 23. (Color online) The multi-boson distorted 4-pion 
correlation function compared to the undistorted case. Also 
shown is the built correlation function constructed from the 
distorted 2-pion correlation function in Fig. 1211 The setting 
with n = 700 roughly correspond to 0 — 5% Pb-Pb collisions 
at the LHC. 


which exploits the expected feature of coherent pion ra¬ 
diation. Pion coherence suppresses QS correlations and 
increases substantially for higher order correlation func¬ 
tions, making multi-pion QS correlations a sensitive mea¬ 
sure of the DCC. An additional possible signature of the 
DCC may lie in the event-by-event 2"'^ order flow vector 
magnitude, q 2 , as defined in Refs. [11,14^. As the DCC 
is expected to radiate near the final state where the ini¬ 
tial spatial anisotropy in non-central heavy-ion collisions 
has greatly diminished, one may expect a different az¬ 
imuthal distribution of coherent pions as compared to 
chaotic pions. In particular, one may expect ~ 0. In 
such a case one expects an anti-correlation between the 
charged pion q 2 and G (coherent fraction). For events 
where the DCC radiates into charged pions, <72 will be 
diluted. Whereas for events of neutral pion radiation, 52 
will be unchanged. 
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The distortions from multi-pion symmetrizations at 
large phase-space densities widens all orders of correla¬ 
tion functions calculated. The distorted built correla¬ 
tion function is also widened to a very similar extent. 
The similarity of the distortions makes the comparison 
of built and measured correlation functions robust to the 
multi-pion effects at moderate densities. 

The early searches for the DCC have not provided any 
clear indication for its existence In those studies, 

the main experimental observable considered was event- 
by-event fluctuations of neutral to charged pion produc¬ 
tion. However, the experimental reconstruction efficiency 
for TT*^ at low pt as compared to charged pions is often 
far too low to allow for precision measurements of such 
fluctuations. As the DCC is expected to radiate at low 
Pt, one also expects an excess in the single-pion spec¬ 
tra at low pt. Measurements of low px pion production 
in heavy-ion collisions have not revealed a dramatic en¬ 
hancement [H, 13 but do not rule out enhancements at 
the 20% level. 

A more promising channel of search for the DCC is one 


Appendix A: 4 -pion plane-wave function 

The fully symmetrized 4-pion plane-wave function is 
given by 

Vpa = ,/J_fe*Gl2^12+'?132:i3-|-923“^23+'J41^41-l-942“^42+9432:43)/4 

V 24^ 

_|_ eh9l2^21+9l3^23 + 'J23^13-l-94i'“42+942^4l4-9434:43)/4 

-|- 5 permutations 

_|_ g»(9l2^21+9l34:24 + 9234:i4-l-94l3'32+942 4:31+943“^34)/4 

-|- 2 permutations 

_|_ g*(9l2^23-l-9l3a:21+923^31-l-94l2'42+942^43+9434^4l)/4 

-|- 7 permutations 

_|_ gh 9 l 2 ^ 41 -l- 9 'i 3 a: 42 + 923 ^ 12 -l- 9 b ^“34 + 942 ^ 31 + 9432:32 )/4 

-|- 5 permutations], (Al) 

where the second (2-pion), third (2-pair), fourth (3- 
pion), and fifth (4-pion) plane-wave represent the differ¬ 
ent symmetrization sequences. 
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